library(foreign)
library(readstata13)
library(stargazer)
library(sjPlot)
library(ggplot2)
library(ggpubr)
library(sampleSelection)
library(mediation)


rm(list=ls())

# Assign the dataset from dta
arabbar2 <- read.dta13("Arab Barometer II_egypt and tunisia v4.dta", convert.factors = FALSE)



# Further data work

arabbar2$supportproteg <- ifelse(arabbar2$eg801==2, 1, 0)
arabbar2$supportproteg[arabbar2$eg801 >3] <- NA

arabbar2$supportprottn <- ifelse(arabbar2$t901==2, 1, 0)
arabbar2$supportprottn[arabbar2$t901 >3] <- NA

arabbar2$supportprot[arabbar2$supportproteg==1 | arabbar2$supportprottn==1] <- 1
arabbar2$supportprot[arabbar2$supportproteg==0 | arabbar2$supportprottn==0] <- 0

table(arabbar2$supportproteg)
table(arabbar2$supportprottn)
table(arabbar2$supportprot)


arabbar2$opportunityprev <- 0- arabbar2$opportunityp

arabbar2$grievances <- arabbar2$econgrie + arabbar2$polsgrie
arabbar2$grievancesegy <- arabbar2$econgrieegy + arabbar2$polsgrie1egy
arabbar2$grievancestun <- arabbar2$econgrietun + arabbar2$polsgrie1tun

arabbar2$educunemp <- ifelse(arabbar2$collegeeduc==1 & arabbar2$employed==0, 1, 0)
arabbar2$uneducunemp <- ifelse(arabbar2$collegeeduc==0 & arabbar2$employed==0, 1, 0)


arabbar2$followgeneral <- (arabbar2$follownews + arabbar2$follownewstv + 
                             arabbar2$follownewspress + arabbar2$follownewsradio + 
                             arabbar2$followinternet) / 5
arabbar2$followgeneral1 <- (arabbar2$follownews + arabbar2$follownewstv + 
                              arabbar2$follownewspress) / 3
arabbar2$followgeneral2 <- (arabbar2$follownewsradio + 
                              arabbar2$followinternet) / 2
arabbar2$followgeneral3 <- (arabbar2$follownewstv + 
                              arabbar2$follownewspress + arabbar2$follownewsradio + 
                              arabbar2$followinternet) / 4
arabbar2$follownewsinternet <- (arabbar2$follownews + 
                                  arabbar2$followinternet) / 2

arabbar2$friendprotesttn <- ifelse(arabbar2$t905==1, 1, 0)
arabbar2$friendprotesttn[arabbar2$t905 >3 | arabbar2$t905==0] <- NA

arabbar2$friendprotestegy <- ifelse(arabbar2$eg806==1, 1, 0)
arabbar2$friendprotestegy[arabbar2$eg806 >3 | arabbar2$eg806==0] <- NA

arabbar2$friendprotest[arabbar2$friendprotesttn==1 | 
                         arabbar2$friendprotestegy==1] <- 1
arabbar2$friendprotest[arabbar2$friendprotesttn==0 | 
                         arabbar2$friendprotestegy==0] <- 0

summary(arabbar2$econgrie)
summary(arabbar2$polsgrie)
summary(arabbar2$grievances)
summary(arabbar2$grievancesegy)
summary(arabbar2$grievancestun)


egarabbar2 <- subset(arabbar2, egypt==1)
tnarabbar2 <- subset(arabbar2, egypt==0)


######################################
######## Models used in Paper ######## 
######################################


# Main Models

log.econpolopp.egy <- glm(protest ~ econgrie + polsgrie + opportunityprc + 
                            age + female + collegeeduc + employed + urban + 
                            interestpols + follownews + membership + 
                            trust + religious + supportdemos +
                            internetuse  
                          , family = binomial(link = "logit") , 
                          weights=wt, data=arabbar2, subset = (egypt=="1"))
summary(log.econpolopp.egy)
log.econpolopp.tun <- glm(protest ~ econgrie + polsgrie + opportunityprc + 
                            age + female + collegeeduc + employed + urban + 
                            interestpols + follownews + membership + 
                            trust + religious + supportdemos +
                            internetuse
                          , family = binomial(link = "logit") , 
                          weights=wt, data=arabbar2, subset = (egypt=="0"))
summary(log.econpolopp.tun)

log.grieopp.egy <- glm(protest ~ grievances + opportunityprc + 
                         age + female + collegeeduc + employed + urban + 
                         interestpols + follownews + membership + 
                         trust + religious + supportdemos +
                         internetuse
                       , family = binomial(link = "logit") , 
                       weights=wt, data=arabbar2, subset = (egypt=="1"))
summary(log.grieopp.egy)

log.grieopp.tun <- glm(protest ~ grievances + opportunityprc + 
                         age + female + collegeeduc + employed + urban + 
                         interestpols + follownews + membership + 
                         trust + religious + supportdemos +
                         internetuse
                       , family = binomial(link = "logit") , 
                       weights=wt, data=arabbar2, subset = (egypt=="0"))
summary(log.grieopp.tun)

log.grieoppint.egy <- glm(protest ~ grievances * opportunityprc + 
                         age + female + collegeeduc + employed + urban + 
                         interestpols + follownews + membership + 
                         trust + religious + supportdemos +
                         internetuse
                       , family = binomial(link = "logit") , 
                       weights=wt, data=arabbar2, subset = (egypt=="1"))
summary(log.grieoppint.egy)

log.grieoppint.tun <- glm(protest ~ grievances * opportunityprc + 
                         age + female + collegeeduc + employed + urban + 
                         interestpols + follownews + membership + 
                         trust + religious + supportdemos +
                         internetuse
                       , family = binomial(link = "logit") , 
                       weights=wt, data=arabbar2, subset = (egypt=="0"))
summary(log.grieoppint.tun)



#### Plotting Grievances and Opportunities (Figure 1)
# Economic Grievances Plot
modeldata_econ_egy <- get_model_data(log.econpolopp.egy,
                                     type = "pred", terms = "econgrie")
modeldata_econ_tun <- get_model_data(log.econpolopp.tun,
                                     type = "pred", terms = "econgrie")


plot_econ_egy <- ggplot(modeldata_econ_egy, aes(x=x))
plot_econ_egy <- plot_econ_egy + geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
                                             alpha=0.1,fill="blue")
plot_econ_egy <- plot_econ_egy + geom_line(aes(x = x, y = predicted), color="blue")
plot_econ_egy <- plot_econ_egy + geom_point(aes(x = x, y = predicted),
                                            lwd = 1, shape = 19, color = "blue")
plot_econ_egy <- plot_econ_egy + xlab("Economic Grievances") + ylab("Predicted Protest")
plot_econ_egy <- plot_econ_egy + scale_y_continuous(limits=c(0,0.53),
                                                    breaks = c(0,0.1,0.2,0.3,0.4,0.5))
plot_econ_egy <- plot_econ_egy + ggtitle("Egypt")
plot_econ_egy

plot_econ_tun <- ggplot(modeldata_econ_tun, aes(x=x))
plot_econ_tun <- plot_econ_tun + geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
                                             alpha=0.1,fill="red")
plot_econ_tun <- plot_econ_tun + geom_line(aes(x = x, y = predicted), color="red")
plot_econ_tun <- plot_econ_tun + geom_point(aes(x = x, y = predicted),
                                            lwd = 1, shape = 19, color = "red")
plot_econ_tun <- plot_econ_tun + xlab("Economic Grievances") + ylab("Predicted Protest")
plot_econ_tun <- plot_econ_tun + scale_y_continuous(limits=c(0,0.53),
                                                    breaks = c(0,0.1,0.2,0.3,0.4,0.5))
plot_econ_tun <- plot_econ_tun + ggtitle("Tunisia")
plot_econ_tun

ggarrange(plot_econ_egy, plot_econ_tun, nrow=1, ncol=2)


# Political Grievances Plot
modeldata_pols_egy <- get_model_data(log.econpolopp.egy,
                                     type = "pred", terms = "polsgrie")
modeldata_pols_tun <- get_model_data(log.econpolopp.tun,
                                     type = "pred", terms = "polsgrie")

plot_pols_egy <- ggplot(modeldata_pols_egy, aes(x=x))
plot_pols_egy <- plot_pols_egy + geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
                                             alpha=0.1,fill="blue")
plot_pols_egy <- plot_pols_egy + geom_line(aes(x = x, y = predicted), color="blue")
plot_pols_egy <- plot_pols_egy + geom_point(aes(x = x, y = predicted),
                                            lwd = 1, shape = 19, color = "blue")
plot_pols_egy <- plot_pols_egy + xlab("Political Grievances") + ylab("Predicted Protest")
plot_pols_egy <- plot_pols_egy + scale_y_continuous(limits=c(0,0.47),
                                                    breaks = c(0,0.1,0.2,0.3,0.4))
plot_pols_egy <- plot_pols_egy + ggtitle("Egypt")
plot_pols_egy

plot_pols_tun <- ggplot(modeldata_pols_tun, aes(x=x))
plot_pols_tun <- plot_pols_tun + geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
                                             alpha=0.1,fill="red")
plot_pols_tun <- plot_pols_tun + geom_line(aes(x = x, y = predicted), color="red")
plot_pols_tun <- plot_pols_tun + geom_point(aes(x = x, y = predicted),
                                            lwd = 1, shape = 19, color = "red")
plot_pols_tun <- plot_pols_tun + xlab("Political Grievances") + ylab("Predicted Protest")
plot_pols_tun <- plot_pols_tun + scale_y_continuous(limits=c(0,0.47),
                                                    breaks = c(0,0.1,0.2,0.3,0.4))
plot_pols_tun <- plot_pols_tun + ggtitle("Tunisia")
plot_pols_tun

ggarrange(plot_pols_egy, plot_pols_tun, nrow=1, ncol=2)


# Combined Grievances Plot
modeldata_grie_egy <- get_model_data(log.grieopp.egy,
                                     type = "pred", terms = "grievances")
modeldata_grie_tun <- get_model_data(log.grieopp.tun,
                                     type = "pred", terms = "grievances")

plot_grie_egy <- ggplot(modeldata_grie_egy, aes(x=x))
plot_grie_egy <- plot_grie_egy + geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
                                             alpha=0.1,fill="blue")
plot_grie_egy <- plot_grie_egy + geom_line(aes(x = x, y = predicted), color="blue")
plot_grie_egy <- plot_grie_egy + geom_point(aes(x = x, y = predicted),
                                            lwd = 1, shape = 19, color = "blue")
plot_grie_egy <- plot_grie_egy + xlab("Economic and Political Grievances") + ylab("Predicted Protest")
plot_grie_egy <- plot_grie_egy + scale_y_continuous(limits=c(0,0.62),
                                                    breaks = c(0,0.1,0.2,0.3,0.4, 0.5, 0.6))
plot_grie_egy <- plot_grie_egy + ggtitle("Egypt")
plot_grie_egy

plot_grie_tun <- ggplot(modeldata_grie_tun, aes(x=x))
plot_grie_tun <- plot_grie_tun + geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
                                             alpha=0.1,fill="red")
plot_grie_tun <- plot_grie_tun + geom_line(aes(x = x, y = predicted), color="red")
plot_grie_tun <- plot_grie_tun + geom_point(aes(x = x, y = predicted),
                                            lwd = 1, shape = 19, color = "red")
plot_grie_tun <- plot_grie_tun + xlab("Economic and Political Grievances") + ylab("Predicted Protest")
plot_grie_tun <- plot_grie_tun + scale_y_continuous(limits=c(0,0.62),
                                                    breaks = c(0,0.1,0.2,0.3,0.4, 0.5, 0.6))
plot_grie_tun <- plot_grie_tun + ggtitle("Tunisia")
plot_grie_tun

ggarrange(plot_grie_egy, plot_grie_tun, nrow=1, ncol=2)


# Opportunity Perception Plot
modeldata_opp_egy <- get_model_data(log.econpolopp.egy,
                                    type = "pred", terms = "opportunityprc")
modeldata_opp_tun <- get_model_data(log.econpolopp.tun,
                                    type = "pred", terms = "opportunityprc")


plot_opp_egy <- ggplot(modeldata_opp_egy, aes(x=x))
plot_opp_egy <- plot_opp_egy + geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
                                           alpha=0.1,fill="blue")
plot_opp_egy <- plot_opp_egy + geom_line(aes(x = x, y = predicted), color="blue")
plot_opp_egy <- plot_opp_egy + geom_point(aes(x = x, y = predicted),
                                          lwd = 1, shape = 19, color = "blue")
plot_opp_egy <- plot_opp_egy + xlab("Opportunity Perception") + ylab("Predicted Protest")
plot_opp_egy <- plot_opp_egy + scale_y_continuous(limits=c(0,0.15),
                                                  breaks = c(0,0.05,0.1,0.15))
plot_opp_egy <- plot_opp_egy + ggtitle("Egypt")
plot_opp_egy

plot_opp_tun <- ggplot(modeldata_opp_tun, aes(x=x))
plot_opp_tun <- plot_opp_tun + geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
                                           alpha=0.1,fill="red")
plot_opp_tun <- plot_opp_tun + geom_line(aes(x = x, y = predicted), color="red")
plot_opp_tun <- plot_opp_tun + geom_point(aes(x = x, y = predicted),
                                          lwd = 1, shape = 19, color = "red")
plot_opp_tun <- plot_opp_tun + xlab("Opportunity Perception") + ylab("Predicted Protest")
plot_opp_tun <- plot_opp_tun + scale_y_continuous(limits=c(0,0.32),
                                                  breaks = c(0,0.1,0.2,0.3))
plot_opp_tun <- plot_opp_tun + ggtitle("Tunisia")
plot_opp_tun


## Figure 1
plot_all_sig <- ggarrange(plot_econ_tun, plot_pols_tun, plot_grie_tun, plot_opp_egy, 
                          nrow=2, ncol=2)
plot_all_sig



#### Regression Table for Main Models (Table 2)

stargazer(log.grieopp.egy, log.grieopp.tun, log.econpolopp.egy, 
          log.econpolopp.tun, log.grieoppint.egy, log.grieoppint.tun,
          style = "apsr",
          covariate.labels = c("Grievances", "Economic Grievances",
                               "Political Grievances", "Opportunity", 
                               "Grievances * Opportunity"),
          dep.var.labels   = c("Protest Participation"),
          object.names = FALSE,
          model.numbers          = TRUE,
          column.labels = c("Egypt", "Tunisia", "Egypt", "Tunisia", "Egypt", "Tunisia"),
          star.char = c("+", "*", "**", "***"),
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
          notes.append = F,
          single.row = FALSE,
          omit = c("age", "female", "collegeeduc", "employed", "urban", 
                   "interestpols", "follownews", "membership", 
                   "trust", "religious", "supportdemos",
                   "internetuse"),
          add.lines = list(c("Controls", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")),
          type="html",
          out="Table2_short.html",
          omit.stat=c("chi2","ser","f"))

stargazer(log.grieopp.egy, log.grieopp.tun, log.econpolopp.egy, 
          log.econpolopp.tun, log.grieoppint.egy, log.grieoppint.tun,
          style = "apsr",
          order = c(1,2,3,4,17,5,6,7,8,9,10,11,12,13,14,15,16),
          covariate.labels = c("Grievances", "Economic Grievances",
                               "Political Grievances", "Opportunity", 
                               "Grievances * Opportunity",
                               "Age", "Female", "College Education", "Employed", "Urban",
                               "Interest in Politics", "Follow News",
                               "Membership", "Trust", "Religiosity", 
                               "Support Democracy", "Internet Use"),
          dep.var.labels   = c("Protest Participation"),
          object.names = FALSE,
          model.numbers          = TRUE,
          column.labels = c("Egypt", "Tunisia", "Egypt", "Tunisia", "Egypt", "Tunisia"),
          star.char = c("+", "*", "**", "***"),
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
          notes.append = F,
          single.row = FALSE,
          type="html",
          out="Table2_full.html",
          omit.stat=c("chi2","ser","f"))



# Opportunity Perception in Egpyt

log.griebycountry <- glm(protest ~ grievances * egypt + opportunityprc +
                           age + female + collegeeduc + employed + 
                           interestpols + follownews + membership + 
                           trust + religious + supportdemos +
                           gradual + urban  + internetuse
                         , family = binomial(link = "logit") , 
                         weights=wt, data=arabbar2)

log.oppbycountry <- glm(protest ~ grievances + opportunityprc * egypt+ 
                          age + female + collegeeduc + employed + 
                          interestpols + follownews + membership + 
                          trust + religious + supportdemos +
                          gradual + urban  + internetuse
                        , family = binomial(link = "logit") , 
                        weights=wt, data=arabbar2)

lm.oppfollow.egy <- lm(opportunityprc ~ follownews + followinternet + follownewstv + follownewspress + follownewsradio +
                         age + female + collegeeduc + employed + urban +
                         interestpols + membership + 
                         trust + religious + supportdemos +
                         gradual + internetuse, 
                       data = arabbar2, weights = wt, subset = (egypt=="1"))
summary(lm.oppfollow.egy)

# Opportunity Perception Plot Egypt /(Figure 2)
modeldata_followopp_egy <- get_model_data(lm.oppfollow.egy,
                                          type = "pred", terms = "follownews")


plot_followopp_egy <- ggplot(modeldata_followopp_egy, aes(x=x))
plot_followopp_egy <- plot_followopp_egy + geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
                                                       alpha=0.1,fill="red")
plot_followopp_egy <- plot_followopp_egy + geom_line(aes(x = x, y = predicted), color="red")
plot_followopp_egy <- plot_followopp_egy + geom_point(aes(x = x, y = predicted),
                                                      lwd = 1, shape = 19, color = "red")
plot_followopp_egy <- plot_followopp_egy + xlab("Follow News") + ylab("Predicted Opportunity Perception")
plot_followopp_egy <- plot_followopp_egy + ggtitle("Egypt")
plot_followopp_egy


# Regression Table Opportunity Perception  (Table 3)

stargazer(log.griebycountry, log.oppbycountry, lm.oppfollow.egy,
          style = "apsr",
          order = c(1,3,2,21,22,9,10,11,12,13,4,5,6,7,14,15,16,17,18,19,20),
          covariate.labels = c("Grievances", "Opportunity Perception", "Egypt",
                               "Grievances * Egypt",
                               "Opportunity Perception * Egypt", "Follow News", "Follow on Internet",
                               "Follow on TV", "Follow on Press", "Follow on Radio"), 
          dep.var.labels   = c("Protest Participation", "Opportunity Perception"),
          object.names = FALSE,
          model.numbers          = TRUE,
          column.labels = c("Both Countries", "Both Countries", "Egypt"),
          star.char = c("+", "*", "**", "***"),
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
          notes.append = F,
          single.row = FALSE,
          omit = c("age", "female", "collegeeduc", "employed", "urban", 
                   "interestpols", "membership", 
                   "trust", "religious", "supportdemos",
                   "internetuse", "gradual"),
          add.lines = list(c("Controls", "Yes", "Yes", "Yes")),
          type="html",
          out="Table3_short.html",
          omit.stat=c("chi2","ser","f", "aic", "adj.rsq"))


stargazer(log.griebycountry, log.oppbycountry, lm.oppfollow.egy,
          style = "apsr",
          order = c(1,3,2,21,22,9,10,11,12,13,4,5,6,7,14,15,16,17,18,19,20),
          covariate.labels = c("Grievances", "Opportunity Perception", "Egypt",
                               "Grievances * Egypt",
                               "Opportunity Perception * Egypt", "Follow News", "Follow on Internet",
                               "Follow on TV", "Follow on Press", "Follow on Radio", 
                               "Age", "Female", "College Education", "Employed", "Membership",
                               "Trust", "Religiosity", 
                               "Support Democracy", "Gradual Change", "Urban",
                               "Internet Use",  "Interest in Politics"), 
          dep.var.labels   = c("Protest Participation", "Opportunity Perception"),
          object.names = FALSE,
          model.numbers          = TRUE,
          column.labels = c("Both Countries", "Both Countries", "Egypt"),
          star.char = c("+", "*", "**", "***"),
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
          notes.append = F,
          single.row = FALSE,
          type="html",
          out="Table3_full.html",
          omit.stat=c("chi2","ser","f", "aic", "adj.rsq"))




########################
###### Appendix ########
######################## 


## Extra Controls (Table A2)

log.control.grieopp.egy <- glm(protest ~ grievances + opportunityprc + 
                                 income + friendprotest + sharia + gradual + 
                                 age + female + collegeeduc + employed + 
                                 interestpols + follownews + membership + 
                                 trust + religious + supportdemos +
                                 urban + internetuse
                               , family = binomial(link = "logit") , 
                               weights=wt, data=arabbar2, subset = (egypt=="1"))
summary(log.control.grieopp.egy)

log.control.grieopp.tn <- glm(protest ~ grievances + opportunityprc + 
                                income + friendprotest + sharia + gradual +                  
                                age + female + collegeeduc + employed + 
                                interestpols + follownews + membership + 
                                trust + religious + supportdemos +
                                urban + internetuse
                              , family = binomial(link = "logit") , 
                              weights=wt, data=arabbar2, subset = (egypt=="0"))
summary(log.control.grieopp.tn)

log.control.econpolsopp.egy <- glm(protest ~ econgrie + polsgrie + opportunityprc + 
                                     income + friendprotest + sharia + gradual + 
                                     age + female + collegeeduc + employed + 
                                     interestpols + follownews + membership + 
                                     trust + religious + supportdemos +
                                     urban + internetuse
                                   , family = binomial(link = "logit") , 
                                   weights=wt, data=arabbar2, subset = (egypt=="1"))
summary(log.control.econpolsopp.egy)

log.control.econpolsopp.tn <- glm(protest ~ econgrie + polsgrie + opportunityprc + 
                                    income + friendprotest + sharia + gradual +                  
                                    age + female + collegeeduc + employed + 
                                    interestpols + follownews + membership + 
                                    trust + religious + supportdemos +
                                    urban + internetuse
                                  , family = binomial(link = "logit") , 
                                  weights=wt, data=arabbar2, subset = (egypt=="0"))
summary(log.control.econpolsopp.tn)



stargazer(log.control.grieopp.egy, log.control.grieopp.tn,
          log.control.econpolsopp.egy, log.control.econpolsopp.tn,
          style = "apsr",
          covariate.labels = c("Grievances", "Economic Grievances",
                               "Political Grievances",  "Opportunity", 
                               "Income", "Friends Protested", 
                               "Sharia Support", "Gradual Change",
                               "Age", "Female", "College Education", "Employed", 
                               "Interest in Politics", "Follow News",
                               "Membership", "Trust", "Religiosity", 
                               "Support Democracy", "Urban", "Internet Use"),
          dep.var.labels   = c("Protest Participation"),
          object.names = FALSE,
          model.numbers          = TRUE,
          column.labels = c("Egypt", "Tunisia", "Egypt", "Tunisia"),
          star.char = c("+", "*", "**", "***"),
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
          notes.append = F,
          single.row = FALSE,
          type="html",
          out="TableA1.html",
          omit.stat=c("chi2","ser","f"))





### Protest Web (Table A3)

log.web.econpolopp.egy <- glm(protestweb ~ econgrie + polsgrie + opportunityprc + 
                                age + female + collegeeduc + employed + 
                                interestpols + follownews + membership + 
                                trust + religious + supportdemos +
                                gradual + urban  
                              , family = binomial(link = "logit") , 
                              weights=wt, data=arabbar2, subset = (egypt=="1"))
summary(log.web.econpolopp.egy)

log.web.econpolopp.tn <- glm(protestweb ~ econgrie + polsgrie + opportunityprc + 
                               age + female + collegeeduc + employed + 
                               interestpols + follownews + membership + 
                               trust + religious + supportdemos +
                               gradual + urban
                             , family = binomial(link = "logit") , 
                             weights=wt, data=arabbar2, subset = (egypt=="0"))
summary(log.web.econpolopp.tn)

log.web.grieopp.egy <- glm(protestweb ~ grievances + opportunityprc + 
                             age + female + collegeeduc + employed + 
                             interestpols + follownews + membership + 
                             trust + religious + supportdemos +
                             gradual + urban
                           , family = binomial(link = "logit") , 
                           weights=wt, data=arabbar2, subset = (egypt=="1"))
summary(log.web.grieopp.egy)

log.web.grieopp.tn <- glm(protestweb ~ grievances + opportunityprc + 
                            age + female + collegeeduc + employed + 
                            interestpols + follownews + membership + 
                            trust + religious + supportdemos +
                            gradual + urban
                          , family = binomial(link = "logit") , 
                          weights=wt, data=arabbar2, subset = (egypt=="0"))
summary(log.web.grieopp.tn)



stargazer(log.web.grieopp.egy, log.web.grieopp.tn, 
          log.web.econpolopp.egy, log.web.econpolopp.tn,
          style = "apsr",
          covariate.labels = c("Grievances", "Economic Grievances",
                               "Political Grievances", "Opportunity", 
                               "Age", "Female", "College Education", "Employed", 
                               "Interest in Politics", "Follow News",
                               "Membership", "Trust", "Religiosity", 
                               "Support Democracy", "Gradual Change",
                               "Urban"),
          dep.var.labels   = c("Online Protest Participation"),
          object.names = FALSE,
          model.numbers          = TRUE,
          column.labels = c("Egypt", "Tunisia", "Egypt", "Tunisia"),
          star.char = c("+", "*", "**", "***"),
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
          notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
          notes.append = F,
          single.row = FALSE,
          type="html",
          out="TableA2.html",
          omit.stat=c("chi2","ser","f"))


#### Plot for All the Models (Figure A1)

plot_all2 <- ggarrange(plot_econ_egy, plot_econ_tun, plot_pols_egy, plot_pols_tun, 
                       plot_grie_egy,  plot_grie_tun, plot_opp_egy, plot_opp_tun,
                       nrow=4, ncol=2)
plot_all2




#### Plot for All Opportunity Perceptions / Follow (Figure A2)

modeldata_followopp_egy <- get_model_data(lm.oppfollow.egy,
                                          type = "pred", terms = "follownews")
modeldata_followintopp_egy <- get_model_data(lm.oppfollow.egy,
                                             type = "pred", terms = "followinternet")
modeldata_followtvopp_egy <- get_model_data(lm.oppfollow.egy,
                                            type = "pred", terms = "follownewstv")
modeldata_followpressopp_egy <- get_model_data(lm.oppfollow.egy,
                                               type = "pred", terms = "follownewspress")
modeldata_followradioopp_egy <- get_model_data(lm.oppfollow.egy,
                                               type = "pred", terms = "follownewsradio")

#
plot_followopp_egy <- ggplot(modeldata_followopp_egy, aes(x=x))
plot_followopp_egy <- plot_followopp_egy + geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
                                                       alpha=0.1,fill="red")
plot_followopp_egy <- plot_followopp_egy + geom_line(aes(x = x, y = predicted), color="red")
plot_followopp_egy <- plot_followopp_egy + geom_point(aes(x = x, y = predicted),
                                                      lwd = 1, shape = 19, color = "red")
plot_followopp_egy <- plot_followopp_egy + xlab("Follow News") + ylab("Predicted Opportunity Perception")
plot_followopp_egy <- plot_followopp_egy + ggtitle("Egypt")
plot_followopp_egy
#
plot_followintopp_egy <- ggplot(modeldata_followintopp_egy, aes(x=x))
plot_followintopp_egy <- plot_followintopp_egy + geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
                                                             alpha=0.1,fill="red")
plot_followintopp_egy <- plot_followintopp_egy + geom_line(aes(x = x, y = predicted), color="red")
plot_followintopp_egy <- plot_followintopp_egy + geom_point(aes(x = x, y = predicted),
                                                            lwd = 1, shape = 19, color = "red")
plot_followintopp_egy <- plot_followintopp_egy + xlab("Follow News (Internet)") + ylab("Predicted Opportunity Perception")
plot_followintopp_egy <- plot_followintopp_egy + ggtitle("Egypt")
plot_followintopp_egy
#
plot_followtvopp_egy <- ggplot(modeldata_followtvopp_egy, aes(x=x))
plot_followtvopp_egy <- plot_followtvopp_egy + geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
                                                           alpha=0.1,fill="red")
plot_followtvopp_egy <- plot_followtvopp_egy + geom_line(aes(x = x, y = predicted), color="red")
plot_followtvopp_egy <- plot_followtvopp_egy + geom_point(aes(x = x, y = predicted),
                                                          lwd = 1, shape = 19, color = "red")
plot_followtvopp_egy <- plot_followtvopp_egy + xlab("Follow News (TV)") + ylab("Predicted Opportunity Perception")
plot_followtvopp_egy <- plot_followtvopp_egy + ggtitle("Egypt")
plot_followtvopp_egy
#
plot_followpressopp_egy <- ggplot(modeldata_followpressopp_egy, aes(x=x))
plot_followpressopp_egy <- plot_followpressopp_egy + geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
                                                                 alpha=0.1,fill="red")
plot_followpressopp_egy <- plot_followpressopp_egy + geom_line(aes(x = x, y = predicted), color="red")
plot_followpressopp_egy <- plot_followpressopp_egy + geom_point(aes(x = x, y = predicted),
                                                                lwd = 1, shape = 19, color = "red")
plot_followpressopp_egy <- plot_followpressopp_egy + xlab("Follow News (Press)") + ylab("Predicted Opportunity Perception")
plot_followpressopp_egy <- plot_followpressopp_egy + ggtitle("Egypt")
plot_followpressopp_egy
#
plot_followradioopp_egy <- ggplot(modeldata_followradioopp_egy, aes(x=x))
plot_followradioopp_egy <- plot_followradioopp_egy + geom_ribbon(aes(x = x, ymin = conf.low, ymax = conf.high), 
                                                                 alpha=0.1,fill="red")
plot_followradioopp_egy <- plot_followradioopp_egy + geom_line(aes(x = x, y = predicted), color="red")
plot_followradioopp_egy <- plot_followradioopp_egy + geom_point(aes(x = x, y = predicted),
                                                                lwd = 1, shape = 19, color = "red")
plot_followradioopp_egy <- plot_followradioopp_egy + xlab("Follow News (Radio)") + ylab("Predicted Opportunity Perception")
plot_followradioopp_egy <- plot_followradioopp_egy + ggtitle("Egypt")
plot_followradioopp_egy


plot_followversions_all <- ggarrange(plot_followintopp_egy,
                                     plot_followtvopp_egy,
                                     plot_followpressopp_egy, 
                                     plot_followradioopp_egy,
                                     nrow=2, ncol=2)
plot_followversions_all

plot_follow_all <- ggarrange(plot_followopp_egy,
                             ggarrange(plot_followintopp_egy,
                                       plot_followtvopp_egy,
                                       plot_followpressopp_egy, 
                                       plot_followradioopp_egy,
                                       nrow=2, ncol=2), 
                             nrow = 2)
plot_follow_all




# Mediation Analysis (Table A5)

med1.egy <-   mediate( lm(opportunityprc ~ grievances +
                            age + female + collegeeduc + employed + urban + 
                            interestpols + follownews + membership + 
                            trust + religious + supportdemos +
                            internetuse, 
                          data=egarabbar2), 
                       glm(protest ~ grievances + opportunityprc +
                             age + female + collegeeduc + employed + urban + 
                             interestpols + follownews + membership + 
                             trust + religious + supportdemos +
                             internetuse, 
                           data=egarabbar2, family = "binomial"), 
                       sims = 1000, boot = TRUE,
                       treat = "grievances", mediator = "opportunityprc", 
                       covariates = list("age" ,"female" , "collegeeduc" , "employed" , 
                                         "urban" , "interestpols" , "follownews" ,
                                         "membership" , "trust" , "religious" , "supportdemos",
                                         "internetuse"),dropobs=T)
summary(med1.egy)


med1.tn <-   mediate( lm(opportunityprc ~ grievances +
                           age + female + collegeeduc + employed + urban + 
                           interestpols + follownews + membership + 
                           trust + religious + supportdemos +
                           internetuse, 
                         data=tnarabbar2), 
                      glm(protest ~ grievances + opportunityprc +
                            age + female + collegeeduc + employed + urban + 
                            interestpols + follownews + membership + 
                            trust + religious + supportdemos +
                            internetuse, 
                          data=tnarabbar2, family = "binomial"), 
                      sims = 1000, boot = TRUE,
                      treat = "grievances", mediator = "opportunityprc", 
                      covariates = list("age" ,"female" , "collegeeduc" , "employed" , 
                                        "urban" , "interestpols" , "follownews" ,
                                        "membership" , "trust" , "religious" , "supportdemos",
                                        "internetuse"),dropobs=T)
summary(med1.tn)


med2.egy <-   mediate( lm(grievances ~ opportunityprc +
                            age + female + collegeeduc + employed + urban + 
                            interestpols + follownews + membership + 
                            trust + religious + supportdemos +
                            internetuse, 
                          data=egarabbar2), 
                       glm(protest ~ grievances + opportunityprc +
                             age + female + collegeeduc + employed + urban + 
                             interestpols + follownews + membership + 
                             trust + religious + supportdemos +
                             internetuse, 
                           data=egarabbar2, family = "binomial"), 
                       sims = 1000, boot = TRUE,
                       treat = "opportunityprc", mediator = "grievances", 
                       covariates = list("age" ,"female" , "collegeeduc" , "employed" , 
                                         "urban" , "interestpols" , "follownews" ,
                                         "membership" , "trust" , "religious" , "supportdemos",
                                         "internetuse"),dropobs=T)
summary(med2.egy)


med2.tn <-   mediate( lm(grievances ~ opportunityprc +
                           age + female + collegeeduc + employed + urban + 
                           interestpols + follownews + membership + 
                           trust + religious + supportdemos +
                           internetuse, 
                         data=tnarabbar2), 
                      glm(protest ~ grievances + opportunityprc +
                            age + female + collegeeduc + employed + urban + 
                            interestpols + follownews + membership + 
                            trust + religious + supportdemos +
                            internetuse, 
                          data=tnarabbar2, family = "binomial"), 
                      sims = 1000, boot = TRUE,
                      treat = "opportunityprc", mediator = "grievances", 
                      covariates = list("age" ,"female" , "collegeeduc" , "employed" , 
                                        "urban" , "interestpols" , "follownews" ,
                                        "membership" , "trust" , "religious" , "supportdemos",
                                        "internetuse"),dropobs=T)
summary(med2.tn)


